Mapping

The motivation for the folliwng excercises is to how to use and show data! How do we create plots to convey information First module is simple maps. ggplot is used well with this software

rm(list = ls())
library(maps)
library(sf)
library(raster)
library(tmap)
library(micromap)
library(ggrepel)
library(ggmap)

Basic Maps

cities <- c('Ashland', 'Corvallis', 'Bend', 'Portland', 'Newport')
longitude <- c(-122.699, -123.275, -121.313, -122.670, -124.054)
latitude <- c(42.189, 44.57, 44.061, 45.523, 44.652)
population <- c(20062, 50297, 61362, 537557, 9603)
dat <- data.frame(cities, longitude, latitude, population)
plot(latitude ~ longitude, data = dat)

#Convert dataframe to simple feature (sf) object The following maps demonstrate that all these maps acan be doen with base r, but there is a better way to do it- ggplot! after this chunk….

dat_sf <- st_as_sf(dat, coords = c("longitude", "latitude"))
st_crs(dat_sf) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
head(dat_sf) #also shows the data
## Simple feature collection with 5 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -124.054 ymin: 42.189 xmax: -121.313 ymax: 45.523
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##      cities population                geometry
## 1   Ashland      20062 POINT (-122.699 42.189)
## 2 Corvallis      50297  POINT (-123.275 44.57)
## 3      Bend      61362 POINT (-121.313 44.061)
## 4  Portland     537557  POINT (-122.67 45.523)
## 5   Newport       9603 POINT (-124.054 44.652)
plot(st_geometry(dat_sf)) 

map('county', region = 'oregon') #maps package is a way to access base maps, good to know
plot(st_geometry(dat_sf), add = T)

#plot points size based on population 
map('county', region = 'oregon')
plot(st_geometry(dat_sf), add = T, pch = 16, cex = sqrt(dat_sf$population * 0.0002))

Now we are going to do it ggplot only the development ggplot has the new geom_sf that is why we need to download off of github Recreating the same plot using geom_sf

#first all spatial files need to be a sf file
state <- st_as_sf(map('county', region = 'oregon', plot = F, fill = T))

library(ggplot2)

ggplot() +
  geom_sf(data = state) +
  geom_sf(data = dat_sf) +
  theme_minimal() #theme_void usually removes everything expect what is plotted but development verison is funky. 

#adding more data!

ggplot(dat_sf) +
  geom_sf(data = state) +
  geom_sf(aes(size = population, colour = population)) +
  theme_minimal() #theme_void usually removes everything expect what is plotted but development verison is funky. 

p <- last_plot() #last plot that was made


#change some of the aestics

p <- p +
  scale_size(range = c(5, 15), guide = F) + #guide just false because developer version legend is bunk. 
  scale_colour_distiller(type = 'div', palette = 9) #this is a way to change colour palette, div is good. 
p

#p <- p + 
#  geom_text(aes(x = longitude, y = latitude, label = cities)) #when making new geoms it is going to look for the data that is specified within the ggplot () 

#best way to deal with points! There is a cool animation that shows how it is done

p <- p + 
  geom_text_repel(aes(x = longitude, y = latitude, label = cities), point.padding = grid::unit(0.5, 'lines'), segment.alpha = 0)

p

Now we are moving to ggmap which utilizes google base maps!

Need to determine bounding box for google map to be happy

Need to mess with zoom value to get the final product in line.

# get the extent
dat_ext <- unname(st_bbox(dat_sf)) #uname strips the names from the boundary mox, st_bbox gets bounding box

# get the base map using the extent
bsmap <-  get_map(location = dat_ext, maptype = 'satellite', zoom = 7)
## Warning: bounding box given to google - spatial extent only approximate.
## converting bounding box to center/zoom specification. (experimental)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=43.856,-122.6835&zoom=7&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
# plot the basemape
ggmap(bsmap)

#this is our starting point to add more data

ggmap(bsmap) +
  geom_sf(data = dat_sf, inherit.aes = F, color = 'red') #use inheriate.aes false to make sure not using the aes from bsmap
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

ggmap(bsmap) +
  geom_sf(data = dat_sf, inherit.aes = F, aes(size = population, colour= population)) +
  geom_text_repel(data = dat_sf, aes(x = longitude, y = latitude, label = cities), point.padding = grid::unit(0.5, 'lines'), segment.alpha = 0, colour = 'white')
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

ggmap(bsmap) +
  geom_sf(data = dat_sf, inherit.aes = F, aes(size = population, colour= population)) +
  geom_text_repel(data = dat_sf, aes(x = longitude, y = latitude, label = cities), point.padding = grid::unit(0.5, 'lines'), segment.alpha = 0, colour = 'white') +
  scale_size(range = c(5, 15), guide = F) + 
  scale_colour_distiller(type = 'div', palette = 9)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Now to do a new excercise

#library(Sf)
library(sp)
library(maps)
library(ggplot2)

data(meuse) #loads data formats that are in packages or are in your directory as an r file
head(meuse)
##        x      y cadmium copper lead zinc  elev       dist   om ffreq soil
## 1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1
## 2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1
## 3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1
## 4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2
## 5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2
## 6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2
##   lime landuse dist.m
## 1    1      Ah     50
## 2    1      Ah     30
## 3    1      Ah    150
## 4    0      Ga    270
## 5    0      Ah    380
## 6    0      Ga    470
meu_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992) #agr = "constant")
#can now use geom_sf! horray!

#ggmap does not work with projected data so... 

meu_sf <- st_transform(meu_sf, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")


ext <- unname(st_bbox(meu_sf))

bsmap_me <-  get_map(location = ext, maptype = 'satellite', zoom = 13)
## Warning: bounding box given to google - spatial extent only approximate.
## converting bounding box to center/zoom specification. (experimental)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=50.974088,5.743115&zoom=13&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
ggmap(bsmap_me)

#this is our starting point to add more data

ggmap(bsmap_me) +
  geom_sf(data = meu_sf, inherit.aes = F, color = 'white') 
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

colnames(meuse)
##  [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"   
##  [8] "dist"    "om"      "ffreq"   "soil"    "lime"    "landuse" "dist.m"
ggmap(bsmap_me) +
  geom_sf(data = meu_sf, inherit.aes = F, aes(color = cadmium))+
    scale_colour_distiller(type = 'div', palette = 9) #type = spectral will work too!
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#yeah!